#include "alignerOutput.h"
#define CHAR_PER_WORD 16

char getChar(int i)
{
	switch(i)
	{
		case(0):
			return 'A';
			break;
		case(1):
			return'C';
			break;
		case(2):
			return'G';
			break;
		case(3):
			return'T';
			break;
	}
}

uint8 getInt(char i)
{
	switch(i)
	{
		case'A':
		case'a':
			return 0;
		case'C':
		case'c':
			return 1;
		case'G':
		case'g':
			return 2;
		case'T':
		case't':
			return 3;
		case(0):
			return 0;
		default:
			printf("unknown base: %c\n", i);
			return 0;	
	}
}

char getComplement(char base) {

	switch(base) {
		case 'a':
			return 't';
		case 'A':
			return 'T';
		case 't':
			return 'a';
		case 'T':
			return 'A';
		case 'c':
			return 'g';
		case 'C':
			return 'G';
		case 'g':
		case 'n':
			return 'c';
		case 'G':
		case 'N':
			return 'C';

		default: {
				 printf("character : %c \n", base);
				 ERROR_EXIT("unknown base");
			 }
	}
}
void reverse_complement(char* outRead, const char* inRead, const uint32 readLength) {
	for(int i = 0; i < readLength; i++) {
		outRead[i] = getComplement(inRead[i]);
	}
	reverse(outRead, outRead + readLength);
}

void getReverseBuffer(char* inBuffer, uint32 readNum, uint32 readLength, char* outBuffer)
{
	int i;
	for(i=0; i<readNum; i++) {
		reverse_complement(outBuffer + i*readLength, inBuffer + i*readLength, readLength);
	}
}

void decompressDNABuffer(uint32* packed, uint32 readNum,uint32 readLength, char* unpacked)
{
	int i,j,k,l;
	int querySizeInWord = (readLength + CHAR_PER_WORD - 1)/ CHAR_PER_WORD; 
	uint32 temp;
	for(i=0; i<readNum; i++) {
		k = 0;
		for(j=0; j<querySizeInWord && k < readLength; j++) {
			temp = packed[i*querySizeInWord + j]; 
			for(l=0; l<CHAR_PER_WORD && k < readLength; l++) {
				unpacked[i * readLength + k] = getChar((temp & 0xC0000000) >> 30);
				temp <<= 2;
				k++;
			}
		}
	}
}
char* getTitle(char* titleBuffer, uint32 readNameLength ,uint32 queryId)
{
	return titleBuffer + readNameLength * queryId;
}

void getQuery(char* queryBuffer, uint32 queryLen, uint32 queryId, char* query)
{
	int i;
	for(i=0; i<queryLen; i++)
		query[i] = queryBuffer[queryId * queryLen + i];
	query[i] = '\0';
}


static void setTail(uint32* size0, uint32 size)
{
	(*size0) = (size / 16 + 1) * 16;
}

void compressInfo(Info* info, uint32 blockSize, FILE* fout)
{
	uint8* h_misNum = info->numMisMatch;

	uint8* c_h_misNum = (uint8*)malloc(ceil((double)blockSize/4.0));
	uint32 compressedSize = 0;
	uint8* h_misPos = (uint8*)malloc(blockSize*2);
	uint8* c_h_misPos = (uint8*)malloc(blockSize * 2);
	uint8* h_varBase = (uint8*)malloc(blockSize*2);
	uint8* c_h_varBase = (uint8*)malloc(blockSize);
	
	uint32 varPosOffset = 0;
	uint32 misPosOffset = 0;
	for(int i=0; i<blockSize; i++) {
		if(h_misNum[i] == 1){
			h_misPos[misPosOffset++] = info->misPos1[i];
			h_varBase[varPosOffset++] = info->misBase1[i];
		}
		else if(h_misNum[i] == 2) {
			h_misPos[misPosOffset++] = info->misPos1[i];
			h_misPos[misPosOffset++] = info->misPos2[i];
			h_varBase[varPosOffset++] = info->misBase1[i];
			h_varBase[varPosOffset++] = info->misBase2[i];
		}
	}

	uint32 c_h_misPosBitSize;
	uint32 c_h_misPosSize, c_varPosOffset;

	uint32 *d_misPos_idx, *d_queryId, *d_queryId_idx; 
	uint8* d_misPos, *d_varBase, *d_misNum,  *c_d_misPos, *c_d_varBase, *c_d_misNum;
	uint32 blockSize0, misPosOffset0, varPosOffset0;
	
	setTail(&blockSize0, blockSize);
	setTail(&misPosOffset0, misPosOffset);
	setTail(&varPosOffset0, varPosOffset);
	
	GPUMALLOC((void**)&d_misPos, misPosOffset0);
	GPUMALLOC((void**)&d_varBase, varPosOffset0);
	GPUMALLOC((void**)&d_misNum, blockSize0);
	
	GPUMALLOC((void**)&c_d_misPos, (misPosOffset0+1) * sizeof(int));
	GPUMALLOC((void**)&d_misPos_idx, (misPosOffset0+1) * sizeof(int));
	GPUMALLOC((void**)&c_d_varBase, varPosOffset0 / 4);
	GPUMALLOC((void**)&c_d_misNum, blockSize0 / 4);

	cutilSafeCall(cudaMemset(d_misPos, 0, misPosOffset0));
	cutilSafeCall(cudaMemset(d_varBase, 0, varPosOffset0));
	cutilSafeCall(cudaMemset(d_misNum, 0, blockSize0));
	cutilSafeCall(cudaMemset(c_d_misPos, 0,misPosOffset0*sizeof(int)));
	cutilSafeCall(cudaMemset(d_misPos_idx, 0,  (misPosOffset0 +1)*sizeof(int)));
	cutilSafeCall(cudaMemset(c_d_varBase, 0, varPosOffset0/4));
	cutilSafeCall(cudaMemset(c_d_misNum, 0, blockSize0/4));
	
	printf("blockSize: %u, misPosOffset: %u, misPosOffset0: %u\n",blockSize, misPosOffset, misPosOffset0);

	TOGPU(d_misNum, h_misNum, blockSize*sizeof(char));
	TOGPU(d_misPos,h_misPos, misPosOffset*sizeof(char));
	TOGPU(d_varBase,h_varBase,varPosOffset *sizeof(char));

	INIT_TIMER;
	START_TIMER;
	uint32 c_misPosOffset = 0;
	uint32 M_misPos = 8;

	gpu_golomb_coding(d_misPos, misPosOffset0, d_misPos_idx, c_d_misPos, &c_h_misPosBitSize, &M_misPos);
	FROMGPU(&c_h_misPosBitSize, d_misPos_idx + misPosOffset, sizeof(int));
	c_h_misPosSize = (c_h_misPosBitSize + 7) / 8;

	printf("misPos before compression: %u, after compression: %u, compression ratio: %.2f\n", misPosOffset, 
			c_h_misPosSize, (double)misPosOffset / (double)c_h_misPosSize);

	gpu_binary_encoding(d_misNum, blockSize0, c_d_misNum);
	//binary compress dna variation
	gpu_binary_encoding(d_varBase, varPosOffset0, c_d_varBase);

	c_varPosOffset = ceil((double) varPosOffset / 4.0);
	
	END_TIMER;
	PRINT_TIMER_SEC("mis info compression on gpu: ");

	FROMGPU(c_h_misPos, c_d_misPos, c_h_misPosSize);
	FROMGPU(c_h_varBase, c_d_varBase, ceil((double)varPosOffset /4.0));
	FROMGPU(c_h_misNum, c_d_misNum, ceil((double)blockSize / 4.0));

	fwrite(&blockSize, sizeof(int), 1, fout);
	fwrite(&M_misPos, sizeof(int), 1, fout);
	fwrite(&misPosOffset, sizeof(int), 1, fout);
	fwrite(&c_h_misPosBitSize, sizeof(uint32), 1, fout);
	fwrite(&varPosOffset, sizeof(int), 1, fout);
	fwrite(&c_varPosOffset, sizeof(int), 1, fout);
	

	fwrite(c_h_misNum, 1, ceil((double)blockSize/4.0), fout);
	fwrite(c_h_misPos,1,c_h_misPosSize, fout);
	fwrite(c_h_varBase,1,c_varPosOffset, fout);

//	END_TIMER
//	PRINT_TIMER_SEC("compress disk access time: ");

	uint32 uncompSize = blockSize * 5;
	compressedSize = ceil((double)blockSize/4.0) + c_h_misPosSize + c_varPosOffset;
	
	printf("mis infomation before compression:%u after compression: %u, compressRatio: %.2f\n",uncompSize, compressedSize, (double)uncompSize/compressedSize);

	free(c_h_misPos);
	free(c_h_varBase);

	free(h_misPos);
	free(h_varBase);

	FREEGPU(d_misNum);
	FREEGPU(c_d_misNum);
	FREEGPU(d_misPos);
	FREEGPU(c_d_misPos);
	FREEGPU(d_misPos_idx);
	FREEGPU(d_varBase);
	FREEGPU(c_d_varBase);

}
void findInfo(Read* readBuf, Buffer* buf, uint32 nresult, uint32 numSA, AlignResult* alignRes, char* ref)
{
	uint32* queryId = buf->queryId;
	uint32* pos = buf->pos;
	uint32* offset = buf->offset;

	uint32 i,j,k;
	uint32* readBuffer = readBuf->readBuffer;
	char* readNameBuffer = readBuf->readNameBuffer;
	int maxReadNameLength = readBuf->maxReadNameLength;
	int readLength = readBuf->readLength;;
	int ReadNum = readBuf->readNum;
	int textLength = readBuf->textLength;
	struct Info info;
	
	info.numMisMatch = (uint8*)malloc(numSA);
	info.misPos1 = (uint8*)malloc(numSA);	
	info.misPos2 = (uint8*)malloc(numSA);	
	info.misBase1 = (uint8*)malloc(numSA);	
	info.misBase2 = (uint8*)malloc(numSA);	

	memset(info.numMisMatch,0, numSA);
	memset(info.misPos1,0, numSA);
	memset(info.misPos2,0, numSA);
	memset(info.misBase1,0, numSA);
	memset(info.misBase2,0, numSA);

	info.readLength = readLength;
	info.numSA = numSA;
	info.readStringLength = readLength * ReadNum * 2;
	info.refStringLength = textLength;

	printf("readNum: %u, readLength: %u, textLength: %u, numSA: %u, numResult: %u\n", ReadNum, readLength, textLength, numSA, nresult);

	char* queryBuffer = (char*)malloc(ReadNum * readLength * 2);	
	decompressDNABuffer(readBuffer, ReadNum*2, readLength, queryBuffer);
//	getReverseBuffer(queryBuffer, ReadNum, readLength, queryBuffer + (readLength * ReadNum));

	findMisMatch(queryBuffer, ref, queryId, pos, offset, &info);	

	alignRes->_numAlign = numSA;
	alignRes->_readLength = readLength;
	alignRes->_pos = pos;

	k=0;
	for(i=0; i<numSA-1; i++) {
		for(j = offset[i]; j<offset[i+1]; j++){
			alignRes->_strand[k] = queryId[i] > ReadNum ? '+' : '-';
			alignRes->_offset[k].x = info.misPos1[k];
			alignRes->_offset[k].y = info.misPos2[k];
			if(info.numMisMatch[k] > 0)
				alignRes->_sub[k].x = info.misBase1[k];
			else 
				alignRes->_sub[k].x = -1;
			if(info.numMisMatch[k] > 1)
				alignRes->_sub[k].y = info.misBase2[k];
			else 
				alignRes->_sub[k].y = -1;
			k++;
		}
	}
	for(; j<nresult; j++){
			alignRes->_strand[k] = queryId[i] > ReadNum ? '+' : '-';
			alignRes->_offset[k].x = info.misPos1[k];
			alignRes->_offset[k].y = info.misPos2[k];
			if(info.numMisMatch[k] > 0)
				alignRes->_sub[k].x = info.misBase1[k];
			else 
				alignRes->_sub[k].x = -1;
			if(info.numMisMatch[k] > 1)
				alignRes->_sub[k].y = info.misBase2[k];
			else 
				alignRes->_sub[k].y = -1;
			k++;
	}

	assert(k == nresult);

	free(info.numMisMatch);
	free(info.misPos1);
	free(info.misPos2);
	free(info.misBase1);
	free(info.misBase2);
	free(queryBuffer);
}



void compressOutputLimited(Read* readBuf, Buffer* buf, uint32 nresult, uint32 numSA, FILE *fout, char* ref)
{
	uint32* queryId = buf->queryId;
	uint32* pos = buf->pos;
	uint32* offset = buf->offset;

	uint32* readBuffer = readBuf->readBuffer;
	char* readNameBuffer = readBuf->readNameBuffer;
	int maxReadNameLength = readBuf->maxReadNameLength;
	int readLength = readBuf->readLength;;
	int ReadNum = readBuf->readNum;
	int textLength = readBuf->textLength;
	struct Info info;
	
	info.numMisMatch = (uint8*)malloc(numSA);
	info.misPos1 = (uint8*)malloc(numSA);	
	info.misPos2 = (uint8*)malloc(numSA);	
	info.misBase1 = (uint8*)malloc(numSA);	
	info.misBase2 = (uint8*)malloc(numSA);	
	memset(info.numMisMatch,0, numSA);
	memset(info.misPos1,0, numSA);
	memset(info.misPos2,0, numSA);
	memset(info.misBase1,0, numSA);
	memset(info.misBase2,0, numSA);
	info.readLength = readLength;
	info.numSA = numSA;
	info.readStringLength = readLength * ReadNum * 2;
	info.refStringLength = textLength;

	printf("readNum: %u, readLength: %u, textLength: %u, numSA: %u, numResult: %u\n", ReadNum, readLength, textLength, numSA, nresult);

	char* queryBuffer = (char*)malloc(ReadNum * readLength * 2);	
	decompressDNABuffer(readBuffer, ReadNum*2, readLength, queryBuffer);
//	getReverseBuffer(queryBuffer, ReadNum, readLength, queryBuffer + (readLength * ReadNum));

/*
	FILE* fout1 = fopen("query", "rb");
	fread(queryBuffer, 1, readLength*2*ReadNum, fout1);
	fclose(fout1);
	
	FILE* fout2 = fopen("pos", "rb");
	fread(pos,sizeof(int),nresult, fout2);
	fclose(fout2);

	FILE* fout3 = fopen("queryId", "rb");
	fread(queryId, sizeof(int), numSA, fout3);
	fclose(fout3);

	FILE* fout4 = fopen("offset", "rb");
	fread(offset, sizeof(int), numSA, fout4);
	fclose(fout4);

*/
	findMisMatch(queryBuffer, ref, queryId, pos, offset, &info);	
		
	compressInfo(&info, numSA, fout);

	free(info.numMisMatch);
	free(info.misPos1);
	free(info.misPos2);
	free(info.misBase1);
	free(info.misBase2);
	free(queryBuffer);
}

void fullSizeOutput(Read* readBuf, Buffer* buf, uint32 nresult, uint32 numSA, FILE *fout, char* ref)
{
	INIT_TIMER;

	uint32* queryId = buf->queryId;
	uint32* pos = buf->pos;
	uint32* offset = buf->offset;

	uint32* readBuffer = readBuf->readBuffer;
	char* readNameBuffer = readBuf->readNameBuffer;
	int maxReadNameLength = readBuf->maxReadNameLength;
	int readLength = readBuf->readLength;;
	int ReadNum = readBuf->readNum;
	int textLength = readBuf->textLength;
	struct Info info;
	
	info.numMisMatch = (uint8*)malloc(numSA);
	info.misPos1 = (uint8*)malloc(numSA);	
	info.misPos2 = (uint8*)malloc(numSA);	
	info.misBase1 = (uint8*)malloc(numSA);	
	info.misBase2 = (uint8*)malloc(numSA);	
	info.readLength = readLength;
	info.numSA = numSA;
	info.readStringLength = readLength * ReadNum * 2;
	info.refStringLength = textLength;

	char* queryBuffer = (char*)malloc(ReadNum * 2 * readLength);
	decompressDNABuffer(readBuffer, ReadNum, readLength, queryBuffer);
	getReverseBuffer(queryBuffer, ReadNum, readLength, queryBuffer + readLength * ReadNum);

	START_TIMER;

	findMisMatch(queryBuffer, ref, queryId, pos, offset, &info);
//	findMisMatch(queryBuffer, ref, queryId, pos, offset, &info);	

	END_TIMER;
	PRINT_TIMER_SEC("findMisMatch: ");

	uint32 i,j,k, id;
	char query[200];
	char *title = ">gnl|SRA|SRR001258.1.1 20563AAXX:8:1:497:903.1 length=47";
	int numMis;
	char strand;
	k = 0;
	START_TIMER;
	fprintf(fout, "title sequence position strand misNum misPos1 misBase1 misPos2 misBase2\n");
	for(i=0; i<numSA-1; i++) {
		id = queryId[i];
		numMis = info.numMisMatch[id];
		strand = '+';
		if(id >= ReadNum) {
			strand = '-';
			id = id - ReadNum;
		}

		getQuery(queryBuffer, readLength, id, query);
//		title = getTitle(readNameBuffer, maxReadNameLength, id);

		for(j=offset[i]; j<offset[i+1]; j++) {
			fprintf(fout, "%s %s %u %c %d ",title,query, pos[k], strand, numMis);
			if(numMis > 0) 
				fprintf(fout,"%d %c ", info.misPos1[queryId[i]], info.misBase1[queryId[i]]);
			if(numMis > 1) 
				fprintf(fout,"%d %c ", info.misPos2[queryId[i]], info.misBase2[queryId[i]]);
			fprintf(fout,"\n");
			k++;
		}
	}
	id = queryId[i];
	numMis = info.numMisMatch[id];
	strand = '+';
	if(id >= ReadNum) {
		strand = '-';
		id = id - ReadNum;
	}

	getQuery(queryBuffer, readLength, id, query);
	title = getTitle(readNameBuffer, maxReadNameLength, id);

	for(j=offset[i]; j<nresult; j++) {
		fprintf(fout, "%s %s %u %c %d ",title,query, pos[k], strand, numMis);
		if(numMis > 0) 
			fprintf(fout,"%d %c ", info.misPos1[queryId[i]], info.misBase1[queryId[i]]);
		if(numMis > 1) 
			fprintf(fout,"%d %c ", info.misPos2[queryId[i]], info.misBase2[queryId[i]]);
		fprintf(fout, "\n");
		k++;
	}
	assert(k == nresult);
	END_TIMER;
	PRINT_TIMER_SEC("write result to disk : ");

	free(info.numMisMatch);
	free(info.misPos1);
	free(info.misPos2);
	free(info.misBase1);
	free(info.misBase2);
	free(queryBuffer);
}	


uint32 getRefPos(uint3* SA_result, uint32 numSA, uint32* saIndex, uint32 numAlign, uint32* queryId,uint32* offset, uint32* pos)
{
	uint32 i,j, k;
	k = 0;
	for(i=0; i<numSA && k <numAlign; i++){
		queryId[i] = SA_result[i].x;
		offset[i] = k;
		for(j=SA_result[i].y; j<=SA_result[i].z && k<numAlign; j++) {
			pos[k] = saIndex[j];
			if(pos[k] == 0)
				printf("k: %d, j: %u, result.y: %u, result.z: %u\n", k, j,SA_result[i].y, SA_result[i].z);
			k++;
		}
	}
	printf("num alignment: %u\n", k);
	return k;
}

void getRefPosResult(uint3* cpuSA, uint32 numSA, uint32 nresult, uint32 textLength, uint32** queryId0,uint32** offset0, uint32** pos0, char* ref_name)
{
	char saFileName[200] = "";// = "/home/luo_group/data/human_ref/ref_good/ref_good.fa.index.sa";
	const char* saFileSuffix = ".sa";
	strcat(saFileName, ref_name);
	strcat(saFileName,saFileSuffix);

	FILE* saFile = fopen(saFileName, "rb");
	if(saFile == NULL) {
		printf("can't open %s\n", saFileName);
	}
	uint32 saInterval;
	uint32 timer = 0, i=0;

	for(i=0; i<6; i++)
		fread(&saInterval, sizeof(int), 1, saFile);

	uint32 saIndexSize = (textLength + saInterval - 1) / saInterval;
	uint32* saIndex = (uint32*)malloc(saIndexSize * sizeof(int));
	uint32* queryId =  (uint32*)malloc(sizeof(int) * numSA);
	uint32* offset =  (uint32*)malloc(sizeof(int) * numSA);
	uint32* pos =  (uint32*)malloc(sizeof(int) * nresult);

	startTimer(&timer);
	fread(saIndex, sizeof(int), saIndexSize, saFile);
	endTimer(&timer, "load sa from disk");
	fclose(saFile);

	startTimer(&timer);
	getRefPos(cpuSA, numSA, saIndex, nresult, queryId, offset, pos);
	endTimer(&timer, "get ref pos ");
	free(saIndex);

	(*queryId0) = queryId;
	(*offset0) = offset;
	(*pos0) = pos;	
}

void compressQueryId(uint32* queryId, uint32 dataSize,uint32 numQuery, uint32* outId, uint32* outLen, uint32* outDataSize)
{
	uint32 *d_queryId, *d_queryId_len, *c_d_queryId, outSize;

	GPUMALLOC((void**)&d_queryId, dataSize * sizeof(int));
	TOGPU(d_queryId, queryId, dataSize * sizeof(int));
	GPUMALLOC((void**)&c_d_queryId, dataSize * sizeof(int));
	GPUMALLOC((void**)&d_queryId_len, dataSize * sizeof(int));

	INIT_TIMER;
	START_TIMER;

	gpu_delta_encoding(d_queryId, dataSize);

	gpu_run_length_encoding(d_queryId, dataSize, c_d_queryId, d_queryId_len, &outSize);

	END_TIMER;
	PRINT_TIMER_SEC("delta + rle kernel: ");

	FROMGPU(outId, c_d_queryId, sizeof(int)*outSize);
	FROMGPU(outLen,d_queryId_len, sizeof(int)*outSize);

	FREEGPU(c_d_queryId);
	FREEGPU(d_queryId_len);
	FREEGPU(d_queryId);

	(*outDataSize) = outSize;
}

void decompressQueryId(uint32* in, uint32* in_len, uint32 dataSize, uint32* out, uint32* outSize)
{
	uint32 maxOutSize = (*outSize);
	uint32* d_out, *d_in, *d_in_len;
	GPUMALLOC((void**)&d_out, maxOutSize * sizeof(int));
	GPUMALLOC((void**)&d_in, dataSize * sizeof(int));
	
	TOGPU(d_in, in, dataSize*sizeof(int));
	GPUMALLOC((void**)&d_in_len, dataSize * sizeof(int));
	TOGPU(d_in_len, in_len, dataSize*sizeof(int));

	INIT_TIMER;
	START_TIMER;

	gpu_run_length_decoding(d_in, d_in_len, dataSize, d_out, outSize);
	gpu_delta_decoding(d_out, *outSize);

	END_TIMER;
	PRINT_TIMER_SEC("delta + rle decoding kernel: ");

	FROMGPU(out, d_out, sizeof(int)*(*outSize));

	FREEGPU(d_in);
	FREEGPU(d_in_len);
	FREEGPU(d_out);
}

void outputIdPos(uint32* c_queryId, uint32* c_queryId_len, uint32 compressedSize, uint32* offset,uint32 numSA, 
		uint32* pos, uint32 nresult, FILE* resultFile, uint32* totalOutputSize)
{
	fwrite(&numSA, sizeof(int), 1, resultFile);
	fwrite(&nresult, sizeof(int), 1, resultFile);
	fwrite(&compressedSize, sizeof(int), 1, resultFile);
	fwrite(c_queryId, sizeof(int), compressedSize, resultFile);
	fwrite(c_queryId_len, sizeof(int), compressedSize, resultFile);
	fwrite(offset, sizeof(int), numSA, resultFile);
	fwrite(pos, sizeof(int), nresult, resultFile); 

	(*totalOutputSize) += ( compressedSize * 2 + numSA + nresult )*4;

	printf("compressed queryid arary: %u bytes\n", compressedSize * 2 * 4);
	printf("offset arary: %u bytes\n", numSA* 4);
	printf("pos arary: %u bytes\n", nresult * 4);
}

void testCompress()
{

	uint32 numSA, nresult,i, timer, textLength;
	const uint32 numQuery = 2000000;
	uint32 compressedSize = numQuery;
	uint32 *c_queryId, *c_queryId_len;
	INIT_TIMER;

	FILE* resultFile = fopen("pos_id_result", "rb");
	if(resultFile == NULL) {
		printf("can't open pos_id_result\n");
		return;
	}

	fread(&nresult, sizeof(int), 1, resultFile);
	fread(&numSA, sizeof(int), 1, resultFile);

	printf("numSA: %u, numAlign: %u\n", numSA, nresult);
	uint32* queryId =  (uint32*)malloc(sizeof(int) * numSA);
	uint32* offset =  (uint32*)malloc(sizeof(int) * numSA);
	uint32* pos =  (uint32*)malloc(sizeof(int) * nresult);

	fread(queryId, sizeof(int), numSA, resultFile);
	fread(offset, sizeof(int), numSA, resultFile);
	fread(pos, sizeof(int), nresult, resultFile);
	/*
	   uint32 k = 0;
	   for(int i =0; i<numSA-1; i++) {
	   for(int j=0; j<(offset[i+1]-offset[i]); j++) {
	   queryId[k++] = queryId_off[i];
	   }
	   }
	   for(int j=0; j<(nresult-offset[numSA-1]); j++) {
	   queryId[k++] = queryId_off[i];
	   }

	   printf("queryId size: %u\n", k);
	   */
	c_queryId = (uint32*)malloc(numQuery * sizeof(int));
	c_queryId_len = (uint32*)malloc(numQuery * sizeof(int));

	START_TIMER;
	compressQueryId(queryId, numSA, numQuery, c_queryId, c_queryId_len, &compressedSize);
	printf("compressed queryId size: %u\n", compressedSize);

	END_TIMER;
	PRINT_TIMER_SEC("rle time for queryId ");

	START_TIMER;
	uint32 *out = (uint32*)malloc(sizeof(int) * nresult);
	uint32 outSize = nresult;
	decompressQueryId(c_queryId, c_queryId_len, compressedSize, out, &outSize);

	END_TIMER;
	PRINT_TIMER_SEC("de-rle time for queryId ");

	for(uint32 i = 0; i<numSA; i++) {
		if(out[i] != queryId[i]) {
			printf("i: %u, out[i]: %u, queryId[i]: %u\n", i, out[i], queryId[i]);
			break;
		}
	}

	free(out);
	free(queryId);
	free(offset);
	free(pos);
}

void testGetRefPos()
{
	uint32 numSA, nresult,i, timer, textLength;
	uint3* cpuSA;
	FILE* saResultFin = fopen("../saInterval", "rb");
	FILE* resultFile = fopen("pos_id_result", "wb");

	char* indexPath = "/home/luo_group/data/human_ref/ref_good/ref_good.fa.index";
	fread(&numSA, sizeof(int), 1, saResultFin);
	fread(&nresult, sizeof(int), 1, saResultFin);
	fread(&textLength, sizeof(int), 1, saResultFin);
	cpuSA = (uint3*)malloc(sizeof(uint3) * numSA);
	fread(cpuSA, sizeof(uint3), numSA, saResultFin);
	fclose(saResultFin);

	printf("numSA: %u, numAlign: %u, textLength: %u\n", numSA, nresult, textLength);

	uint32* queryId, *offset, *pos;
	getRefPosResult(cpuSA, numSA, nresult, textLength, &queryId, &offset, &pos, indexPath);

	printf("write infomation out...\n");
	startTimer(&timer);
	fwrite(&nresult, sizeof(int), 1, resultFile);
	fwrite(&numSA, sizeof(int), 1, resultFile);
	fwrite(queryId, sizeof(int), numSA, resultFile);
	fwrite(offset, sizeof(int), numSA, resultFile);
	fwrite(pos, sizeof(int), nresult, resultFile);
	endTimer(&timer, "write queryId and pos to disk ");

	free(cpuSA);
	free(pos);
	free(offset);
	free(queryId);
	fclose(resultFile);


}

/*
int main()
{

	testCompress();

	return 0;
}

*/
